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The recovery of long-term climate proxy records with seasonal resolution is rare 
because of natural smoothing processes, discontinuities and limitations in 
measurement resolution. Yet insolation forcing, a primary driver of multimillennial- 
scale climate change, acts through seasonal variations with direct impacts on 
seasonal climate’. Whether the sensitivity of seasonal climate to insolation matches 
theoretical predictions has not been assessed over long timescales. Here, we analyse a 
continuous record of water-isotope ratios from the West Antarctic Ice Sheet Divide ice 
coreto reveal summer and winter temperature changes through the last 11,000 years. 


Summer temperatures in West Antarctica increased through the early-to-mid- 
Holocene, reached a peak 4,100 years ago and then decreased to the present. 
Climate model simulations show that these variations primarily reflect changes in 
maximum summer insolation, confirming the general connection between seasonal 
insolation and warming and demonstrating the importance of insolation intensity 
rather than seasonally integrated insolation or season duration??. Winter temperatures 
varied less overall, consistent with predictions from insolation forcing, but also 
fluctuated in the early Holocene, probably owing to changes in meridional heat 
transport. The magnitudes of summer and winter temperature changes constrain the 
lowering ofthe West Antarctic Ice Sheet surface since the early Holocene to less than 
162 m and probably less than 58 m, consistent with geological constraints elsewhere 


in West Antarctica*". 


Milankovitch famously postulated that variations of Earth's orbit and 
axis drive climate changes over tens of thousands of years by altering 
the seasonal cycle of insolation’. By controlling summer temperatures 
and ice ablation, summer insolation in the northern high latitudes is 
thought to drive global ice volume changes over glacial-interglacial 
timescales’. Although modelling studies support this idea?!?, empirical 
evidence of the specific climate response to insolation changes derives 
almost entirely from mean annual temperature reconstructions!"? or 
from indirect effects on, for example, trapped gases and melt layers in 
polar ice?“ and marine aeolian deposits». The absence of seasonal tem- 
perature reconstructions has precluded direct evidence of insolation 
forcing on seasonal climate, a relationship that may vary geographi- 
cally. In Antarctica, long records of multipleglacial-interglacial cycles 
have supported different claims about whether the effects of summer 
insolation relate most strongly to its maximum intensity, its seasonal 
integral or to duration above a threshold??!€7 Site-specific empirical 
determinations would provide valuable tests of such competing ideas. 


Seasonal temperature reconstructions 


Wereconstructed seasonal temperature variability in West Antarctica 
through the Holocene (the last 11,000 years) and performed model 
experiments to understand its physical controls. The Holocene 
offers a window of time for assessing the influence of orbital forcing 
without the complicating effects of Northern Hemisphere deglacia- 
tion”. Our reconstruction (Figs. 1 and 2) uses the high-resolution 
water-isotope record (6D) from the West Antarctic Ice Sheet (WAIS) 
Divide ice core (WDC) (Methods—Water isotopes; Extended Data 
Fig. 1a,b), obtained with a continuous-flow technique that provides 
millimetre-scale depth resolution”. Layer ages were determined pre- 
viously”. 

Records of seasonal temperatures from ice cores are limited by 
measurement resolution and information loss from water-isotope 
diffusion. In Greenland, the longest records separating summer and 
winter variability extend to only 2 thousand years ago (ka) (refs. ?*?*), 
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Fig. 1| Water-isotope seasonal variability. a, Example section of the diffusion- 
corrected (solid line) and raw?? (dashed line) WDC 8D records, with annual 
maxima (redcircles) and minima (blue circles) determined algorithmically 
(Methods-Seasonal water-isotope amplitudes). Extended Data Fig. 1 provides 
the full high-resolution WDC 8D record, diffusion lengths and extrema. 


whereas only climate model simulations are available for older 
periods", For Antarctica, before the present study, the longest records 
spanned only a few centuries”. A combination of three factors accounts 
for the considerably greater scope of our reconstruction: excep- 
tional depth resolution of measurements, conditions at WAIS Divide 
(high accumulation, low temperature and thick ice) which allow for 
preservation of subannual information through the entire Holocene” 
andananalysis strategy which circumvents interannual noise by evalu- 
ating millennial averages ofthe seasonal parameters. 

Our method corrects water-isotope variations for diffusion?*? and 
assesses uncertainties including preservation bias and precipitation 
intermittency (Methods- Diffusion corrections and Uncertainties in 
reconstructing temperatures). The diffusion correction operates on 
the high-resolution data and produces isotopic time series from which 
seasonal summer-winter amplitudes were extracted. These were con- 
verted to temperature using a model-derived scaling? (6.96% 6D C *; 
Methods-Seasonal temperatures) and added to previously recon- 
structed annual mean temperatures? to obtain summer and winter 
histories. 


Age (ka) 


b, The 50-yr annual-amplitude averages (summer minus winter divided by 2), 
with 2c uncertainty; horizontal line indicates Holocene mean. c-e, The 50-yr 
6D averages for summer (red, c), mean (purple, d) and winter (blue, e); 
horizontal line indicates Holocene mean; shaded regions are 20 bounds for 
combined analytical and diffusion-correction uncertainty. 


Seasonal trends 


Summer temperatures at WAIS Divide (Fig. 2a) generally rose through 
the early and middle Holocene, persisted at amaximum between about 
5and 1.5 ka, then decreased toward the present, with a total Holocene 
range of around 2 °C. These variations broadly correlate with local 
maximum insolation, rather than with integrated summer insolation or 
the duration of summer (Fig. 3d,e). Winter temperatures (Fig. 2c) varied 
less than summer ones overall (about 1 °C range) but also fluctuated 
at about 10 to 8 ka, a variation too rapid to attribute to orbital forcing. 

Annual mean WAIS Divide temperature changes” (Fig. 2e) were 
considerably influenced by winter variability in the early Holocene, 
whereas summer variability dominates the overall Holocene pattern 
(Methods-Relationship between the annual mean and individual sea- 
sons; Extended Data Table 2). Summer variability also accounts for 
most of the cooling in the last 2 kyr, indicating that the approximately 
1°C annual-average cooling of the entire West Antarctic during this 
period®” likewise reflects this season. Neither season at WDC experi- 
enced the early Holocene optimum nor overall Holocene cooling that 
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Fig.2|Seasonaltemperature reconstruction. a,c, Reconstructed summer and 
winter temperatures at WDC for 1,000-yr averaging (solid red and blue lines). 
Shaded regions are loand 2ouncertainty ranges for combined uncertainties 
arising from analysis, diffusion correction, seasonality of accumulation, 
precipitation intermittency, isotope-temperature scaling and reconstructed 
mean temperatures (Methods-Uncertainties in reconstructing temperatures). 
Also shown are MEBM-calculated temperatures for 80° S (maximum and 
minimum annual values) and HadCM3 zonal temperatures for 80? S 
(late-December for summer, mid-August for winter) (ORBIT, GLACID and ICE-6C). 
TheO ka ORBIT simulation uses pre-industrial settings, a calculation not 
available for GLACID or ICE-6G. Normalization is done at 1 ka when all model 
runs intersect within 0.05 °C and the ice-sheet configuration is well known. 


appears in some global temperature reconstructions”. To assess the 
significance of the dominant multimillennial trends in each season, 
we performed Monte Carlo analysis (Methods—Trend analysis) using 
4 kaasademarcation pointin summer (thisisthetiming of maximum 
summer temperature) and 6 ka in winter (when winter temperatures 
plateau). For summer (Fig. 2b) this indicates a >95% chance that warm- 
ing from 11to 4 ka and cooling from 4 ka to present exceeded 0.7 and 
0.6 °C, respectively. For winter, the trend from 11 to 6 kais indistinguish- 
able from zero, whereas cooling of greater than about 0.3 °C from 6 to 
0 ka occurred with >95% likelihood (Fig. 2d). 


Moist energy balance model 


To evaluate how orbitally driven insolation changes may explain the 
WAIS Divide reconstructed temperatures (Fig. 2), we first simulated 
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The ICE-6G values at 11 ka for summer and winter (not shown on plots) are 

-3.93 °C and -10.82 °C, respectively. Coefficient of determinations for model 
results versus WDC temperatures (Extended Data Fig. 7) are high for summer 
(HadCM3 ORBIT R?= 0.93, P « 0.001; MEBM A? = 0.80, P « 0.001) but not for 
winter (HadCM3 ORBIT R? = 0.00, P= 0.85; MEBM R? = 0.05, P= 0.30). The winter 
agreement improves if only the period 0-6 ka is considered (HadCM3 ORBIT 

R? =0.74, P=0.01; MEBM R? = 0.39, P= 0.02). b,d, Histograms of net temperature 
changes over the specified time intervals, derived by Monte Carlo analysis 
accounting for systematic and non-systematic uncertainties (Methods—Trend 
analysis).e, WDC mean annual temperature with loand 2ouncertainty bounds”. 
Extended Data Table 2 shows the amount of variability in the mean annual 
temperature that can be explained by the summer and winter temperatures. 


temperature history at 80° S using a global, zonal mean (2° resolution) 
moist energy balance model (MEBM) accounting for incoming and 
outgoing radiation, albedo and meridional atmospheric-heat trans- 
port (Methods—Moist energy balance model). The model is driven 
by top-of-atmosphere (TOA) seasonal insolation changes (Fig. 3a—e); 
for this latitude, the maximum summer insolation increases until 
about 2.5 ka and annual mean and annual- and summer-integrated 
values mostly decline through the Holocene. The calculations 
yield summer maximum temperatures and seasonal temperature 
amplitudes (Fig. 3g) that covary with local maximum summer inso- 
lation (Fig. 3e) and with the general pattern of our reconstructed 
summer temperatures (Extended Data Fig. 7). Although heating at 
lower latitudes can influence Antarctic temperature through atmos- 
pheric and oceanic heat transport, modelled maximum summer 
temperatures at WAIS Divide correlate best with local insolation 


b c 
550 600 600 
_ 500 550 
T D T 
p 9e E 400 E 500 
E = z 
= wae d E = E 
S 500 os 5 300 Z 5 450 É 
E ae E * $ * 
o a 9 200 8 5 400 8 
2 c p] 
£475 — December = E: £ E 
— January 100 2 350 2 
450 Eis ies 0 < 300 < 
11109 8765 43210 123 45 67 8 9 1011 12 11 12 1 2 
Age (ka) Month Month 
d T f 
184 0.10 €E 0 
— 
183 ne a me 
ra] 
j5 0.05 F $ T -20 
= 182 o = E -30 
= 0 S VIE = 
= A = -40 
mt E 
5181 5 9 £5 
3 -005 $ 5 S 5 B» 
2180 9 3535 ce 9 -60 
E & c ` 2 c 
—— Annual mean -010 5 ^ 539l —— Annual maximum 1205 t£ -70 
179} —— Summer (2250 W m^?) "8 — — »275Wm^? Sa 118 5 -20F 
— — Annual © 525+ __ 3950 Wm ~~ Ss] a — —— WDC site -80 
D 
178 » -0.15 © ; 116 -25 -90 
11098765432 E 11109 876543210 11109 87 654 32 10 
Age (ka) Age (ka) Age (ka) 
g MEBM 


Temperature (°C) 
oO 


-0.4 — Maximum summer 
—— Annual amplitude 
—— Minimum winter 
-0.6 
L L L L L L L L L L [ 
11 10 9 8 7 6 5 4 3 2 1 0 
Age (ka) 


Fig.3| Temporal and spatial variability in insolation and model 
temperatures.a, Insolation changethrough the Holocene? for December 
andJanuary and their average. December best resembles the WDC summer 
reconstruction. b, The full seasonal cycle of insolation at 80? S for 500-yr 
snapshots over the Holocene. Line colours inb and c correspond to age. 

c, Zoom of summer insolation. The maximum always occurs in the latter half 
of December (grey shading), migrating across 8 days over the course of the 
Holocene. d, Holocene trends of annual mean insolation (black), annual 


(70° to 90° S, R? = 0.9, P « 0.001 during 0 to 6 ka) rather thaninsolation 
anywhere inthe subtropical through subpolar latitudes (20? to 60? S, 
R? = 0.33-0.55, P< 0.05). Indeed, models indicate heat export from 
WAIS Divide in summer (Extended Data Fig. 4k), rather than import 
from more northern locations. Since December is always the month 
of maximum insolation (Fig. 3a-c), variability of December insola- 
tion dominates the response of maximum summer temperature. For 
winter, modelled temperatures are less variable than those of summer 
at 80? S (Fig. 3g) because of the lack of direct insolation (Fig. 3b) and 
have an opposite trend. Winter minima are a function of three fac- 
tors: changes in the length of the zero-insolation season, the effective 
cooling rate of the surface and convergent heat transport from lower 
latitudes. Lower minimum winter temperatures occur at times when 
the zero-insolation season is longer. However, neither the length 
of the zero-insolation season, modelled minimum temperatures, 
nor winter heat divergence correlate well with reconstructed winter 
temperatures. 


integrated insolation (dashed red line) and summer-integrated insolation 

(red line).e, Maximum summer insolation intensity (black line) and summer 
duration (red lines), defined as the number of days above a threshold insolation 
value each year. f, Anomaly in maximum insolation coloured by latitude in the 
Southern Hemisphere. The thick blue line shows the latitude of the WDC site. 
g, Calculated temperatures for 80° S using the MEBM, including maximum 
summer value (red), minimum winter value (blue) and amplitude of the 
seasonal temperature cycle (black). 


HadCM53 simulations 


Toinvestigate the role of more-complex geography and mechanisms, 
including topographical changes not accounted for in the MEBM, we 
simulated Holocene climate with a fully coupled general circulation 
model, HadCM3 (ref. °°) (Methods—HadCM3 model simulations). Simu- 
lationsforced solely by changes in orbital parameters produce summer 
maximum temperatures (for approximately the December solstice) at 
80° Ssimilartoour reconstructed values and to the MEBM: increasing 
over the Holocene, peaking at 4 to 3 ka and decreasing into the modern 
era (ORBIT, Fig. 2a). This pattern reflects a strong role of maximum 
summer insolation in determining observed summer temperatures. 
The similarity of the early- to mid-Holocene (11-6 ka) summer tem- 
perature increase in the orbitally forced HadCM3 simulations and our 
reconstruction suggests little influence of changing ice-sheet elevation 
and extent. A similar comparison for winter yields an approximately 
1.25 °C decrease of model ORBIT temperature (Fig. 2c) compared to a 
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Fig. 4| Possible ice elevation histories and the corresponding modelled 
temperatures. a, Elevation histories used in HadCM3. GLACID? is 96 m higher 
at 11 ka compared to present and ICE-6G^ is 222 m higher. b, Temperature 
anomalies from elevation change (GLACID, solid lines; ICE-6G, dashed lines) 
using an atmospheric lapse rate of 9.8 °C km ' and spatial lapse rates for 
interior West Antarctica of 12 °C km ! (ref. 5) and 14 °C km7 (ref. ^5). c, HadCM3 


possible small increase in temperature in the reconstruction (Fig. 2d; 
>90% chance of >0.1 °C), suggesting some warming resulting from a 
lowering ice sheet. 

Next, as boundary conditions in the HadCM3 simulations, we pre- 
scribed variable greenhouse gas (GHG) concentrations and two differ- 
ent ice-sheet histories, GLACID and ICE-6G, which entail net surface 
lowerings of about 83 m and about 208 m, respectively, from 11to 7 ka 
at the WDCssite (Fig. 4a). These elevation scenarios substantially affect 
simulated temperatures (Fig. 2a,c). Much of the elevation-induced 
warming in these models, which occurs primarily inthe early Holocene, 
can be attributed directly to the surface lapse-rate effect (Fig. 4b). 
However, comparison to the orbital-only runs (Fig. 4c) reveals a remain- 
ing temperature anomaly (Fig. 4d), attributable to GHGs, ice-sheet 
extent and nonlinear responses to simultaneously imposed forcings. 
Sea ice has only a small impact onthe temperature at 80° Sin summer 
(Methods-Sea ice; Extended Data Fig. 6). 

Inconsistencies exist between the different ice-sheet scenarios 
(Fig. 4d) and the summer versus winter seasons but differences are 
minor enough to permit a bounded estimate of the true Holocene eleva- 
tion decrease. This calculation is made by comparing the excess of the 
reconstructed temperature increase over the orbital-only simulation to 
thesame excess forthe ice-sheet model simulations and scalingto the 
elevation changes used in the latter (Methods—Estimating elevation 
changes). We find central estimates for elevation decrease of 23 m and 
53m from comparison to the GLACID and ICE-6G scenarios, respec- 
tively, over the period 10 to3.5 ka (Table 1). Accounting for uncertainties 
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residual-temperature anomalies for December (summer) calculated by 
subtracting the ORBIT run from the GLACID and ICE-6G runs in Fig. 2a, 
highlighting the portion attributable to changing elevation rather than 
insolation. d, Residual-temperature change in b subtracted from the results 
inc, showing the component driven from processes besides the direct 
lapse-rate (LR) effect and orbital forcing. 


in the seasonal temperature reconstructions (Fig. 2) allows for eleva- 
tion changes ranging from 33 m increase to 131 m decrease (20) from 
10 to 3.5 ka or 54 m increase to 162 m decrease (20) if the time interval 
is narrowed to 10 to 6.5 ka (Table 1). Our results, thus, are consistent 
with geological observations of ice high-stands on mountain nunataks, 
which indicate less than 100 m of Holocene surface lowering * 5. 

Winter temperatures on the Antarctic mainland must respond to 
insolation forcings indirectly, via heat transport from lower latitudes. 
Orbital forcing models predict winter cooling across the Holocene, 
mostly from 11to 6 ka (Figs. 2c and 3g). Both models and reconstructed 
winter temperatures lack alate Holocene maximum. But in the earlier 
Holocene, the winter reconstruction does not display the cooling trend 
expected from models and is dominated by prominent millennial vari- 
ations. The mismatch with insolation at lower latitudes and absence 
of local forcings suggests variations in the efficacy of meridional 
atmospheric-heat transport. 


Discussion 


Diverse and numerous proxies are used to reconstruct globally aver- 
aged surface temperatures for evaluating climate models and distin- 
guishing natural from anthropogenic climate variability?9?*295, How 
these proxies depend on seasonal factors has been assessed in only a 
few cases”. Our West Antarctic study provides a cautionary example, 
asthe mean annual temperature history reflects different control- 
ling factors of summer and winter temperatures whose importance 


Table 1| WAIS elevation decrease 


GLAC1D ICE-6G 
Interval -20 Nominal 20 -20 Nominal 20 
10-6.5ka -9.93 25.67 58.75 -54.00 57.52 161.96 
10-3.5ka -5.63 22.99 48.49 -33.23 52.63 130.53 


Elevation decrease estimates in metres (20, positive values correspond to a lowering ice 
sheet) for the intervals 10-6.5 ka and 10-3.5 ka (Methods—Estimating elevation changes). 


varies with time. In such a situation, important seasonal dynamics 
may be missed, or proxies misinterpreted, when only mean climate 
is considered. In addition, incorporating more information from the 
southern polar regions should help global temperature assessments 
avoid biases associated with weighting oftemperature reconstructions 
toward northernsites, which have produced differing interpretations of 
the relationship between global climate and forcings in the Holocene, 
even including opposing trends?*^?^!, 

Previous analyses with simplified atmospheric models? identified the 
duration of Southern Hemisphere summer as a key driving variable of 
Antarctic climate at orbital timescales. Some palaeoclimate findings 
validate this claim; for example, the onset of deglacial warming in West 
Antarctica corresponds with increasing integrated summer insolation’. 
Our results—spanning about half a precession cycle—reveal a dominant 
role for annual maximum insolation in determining West Antarctic 
summer climate during the Holocene, without precluding a greater 
role for duration or integrated summer insolation in other periods, 
such as glacial terminations. 
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Methods 


We measured WAIS Divide core (WDC) water isotopes using 
continuous-flow analysis (see next section on Water isotopes) and 
then corrected for cumulative diffusion using spectral techniques 
to determine diffusion lengths and restore prediffused ampli- 
tudes within 140-yr sliding windows (section below on Diffusion 
corrections; Extended Data Fig. 1c). Summer maxima and winter 
minima (Fig. 1a) identified in these corrected data were then used 
to calculate summer and winter amplitudes for each year. We con- 
verted the isotope amplitudes to temperature amplitudes using a 
model-determined scaling factor (section below on Seasonal tem- 
peratures) and added them to previously reconstructed mean annual 
temperatures? to recover summer and winter values. Substantial 
seasonal noise processes required multicentennial to millennial 
averaging to reduce uncertainty (section below on Uncertainties 
in reconstructing temperatures). To elucidate physical controls on 
subannual temperatures, we used a simple energy balance model and 
HadCM3, a general circulation model, to calculate expected changes 
in seasonal and monthly surface temperatures through time under 
varying boundary conditions (sections below on Seasonal moist 
energy balance model and HadCM3 model simulations). Finally, using 
both observations and modelling, we estimated the change in WAIS 
surface elevation through the Holocene (section below on Estimating 
elevation changes). 


Water isotopes 

WDC water isotopes (Extended Data Fig. la) were analysed on a 
continuous-flow analysis system” using a Picarro cavity ring-down 
spectroscopy instrument, model L2130-i. Using permutation entropy", 
we identified data anomalies arising from laboratory analysis, which 
were corrected, including by resampling ice through 1,035.4-1,368.2 m 
depths (4,517-6,451 yr)*®. All other Holocene data are previously pub- 
lished'8 and available online???, Data are reported in 5-mm incre- 
ments in delta-notation (%o, or per mil) relative to Vienna standard mean 
ocean water (6/50 = 6D = 096), normalized to standard light Antarctic 
precipitation (6/50 = —55.5%o, 6D = —428.096)). WDC is annually dated, 
with accuracy better than 0.5% of the age between O and 12 ka (ref. ”). 
FortheHolocene, the temporal spacing of consecutive 5 mm samples 
is «0.1 yr and the average «0.05 yr, ranging from about 2.6 weeks at 
10 kato 0.5 week at 1to O ka (ref. !5). 


Diffusion corrections 


Diffusion in the firn and deeper ice attenuates high-frequency water- 
isotope information in ice cores???" ©. Diffusion length quantifies 
the statistical vertical displacement of water molecules from their 
original position”””’. We used diffusion-correction code developed by 
S. Johnsen, University of Copenhagen???*?25. which uses maximum 
entropy methods to invert an observed power-density spectrum. 
As an input to these inversions, we determined diffusion lengths 
(Extended Data Fig. 1c) for 140-yr windows using previous methods'*?6, 
The power-density spectrum observed in the ice-core record P (f), 
after diffusion, is P (f) = P, (f) exp[- Qrfo;) ?j in which P, (f) represents 
the power spectrum of the undiffused signal (96? m”), f is the frequency 
1(/m), Athe signal wavelength (m), zthe depth (m) ando, the diffusion 
length (m). Theoriginal, prediffusion power-density spectrum (diffusion- 
corrected) is calculated as P, (f) = P (f)exp anf? 02) for diffusion length 


0, (yr) and f now with units of I/yr. Theg,= <2 ,in which Awe isthe 
avg 


mean annual layer thickness (m yr“) at a given depth. The diffusion- 
corrected spectrum takes the form of a series of complex numbers 
X + iX, versus f. From this, the amplitude spectrum A is obtained 
by A(f) - X2* X2 and the phase spectrum @ is obtained by 
(f) =tan? B . The real components of the amplitude and phase 
spectrums give the diffusion-corrected water-isotope signal 6, (t) as: 


N 
6,(t) = } A;cos(2nf t+ h) 
i=l 


Uncertainties on 6,(¢) are determined using the uncertainty range 
for diffusion lengths” calculated in each 140-yr window. Before spec- 
tral analysis, the isotope data are linearly interpolated at a uniform 
time interval of 0.05 yr. Our determination of diffusive attenuation 
andcorrectionarises fromthe observed frequency spectra themselves 
and therefore is entirely independent of firn diffusion and densification 
models. 


Seasonal water-isotope amplitudes. To select extrema (summers 
and winters) in the diffusion-corrected 6D signal (Fig. laand Extended 
Data Fig. 1b), we used the ‘findpeaks’ MATLAB function. Figure Ic,e 
show the resulting time series for summer and winter, averaged with 
a50-yr boxcar filter for clarity of trends. For every year defined in the 
WDC age'scale, we calculated the averaged diffusion-corrected 6D. The 
difference between the two extrema and the mean define the summer 
and winter isotope amplitudes. 


Seasonal temperatures. A linear scaling converted seasonal iso- 
topic amplitudes to seasonal temperature amplitudes, using a sensi- 
tivity of isotopes to surface temperatures determined by the simple 
water-isotope model (SWIM)”. Finally, to find summer and winter tem- 
peratures we added the individual seasonal temperature amplitudes 
tothe year's mean temperature obtained previously?? by calibrating 
the water-isotope record against borehole temperatures and 65N con- 
straints on firn thickness. 

SWIM is based on earlier numerical Rayleigh-type distillation 
models****, which simulate the transport and distillation of moisture 
down climatological temperature gradients. As moist air is trans- 
ported towards the poles and cools, the saturated vapour pressure 
decreases nonlinearly and moisture above saturation is removed by 
precipitation. The model keeps track of the isotopic fractionations 
ateach step along this distillation process. In most previous simple 
models, there is an inconsistency in the calculation of the supersat- 
uration that determines the point of condensation and that drives 
kinetic isotope fractionation. Modifications to these earlier models, 
used in SWIM, ensure consistency in the calculation, which results in 
asmoother relationship between temperature and the 6-values of 
precipitation and better agreement with observed spatial patterns of 
6D and 6*0. Given input of both SD and 650 data, SWIM calculates 
distributions of source temperatures, the temperature gradients of 
pseudo-adiabatic pathways and condensation temperature. We used 
SWIM to derive sensitivities for surface isotope-temperature scal- 
ings using diffusion-corrected WDC data to obtain a surface scaling 
of 6.96% 6D °C™. Using raw data, the surface scaling is 7.0726» 8D °C". 
In comparison to other isotope-temperature scalings, ref. ? obtain 
about 6.5696, 6D °C% and ref. °° about 7.1096 6D °C! (both converted 
from 60 to SD using a factor of 8). 


Uncertainties in reconstructing temperatures 

We included uncertainties associated with the following factors: meas- 
urement analysis, diffusion correction, seasonality of accumulation, 
precipitation intermittency, modelled isotope-temperature scaling 
and mean-temperature history. The ‘analysis uncertainty’ is 0.5596» for 
6D (10) (ref. ”). The ‘diffusion-correction uncertainty’ is described in 
ref.?6, The uncertainty of the mean-temperature reconstruction, calcu- 
lated previously??, accounts for mostuncertainty in the early Holocene 
but asmall fraction in the late Holocene. Sections ‘Seasonal preserva- 
tion bias uncertainty’ to 'Isotope-temperature scaling and associated 
uncertainty’ below explain the other uncertainty terms. Uncertainties 
for some factors (analysis and diffusion correction) can be treated as 


independent random variables so that, on time-averaging, their magni- 
tudes decrease as the inverse ofthe square root ofthe number of values. 
Uncertainties for other factors (intermittency, isotope-temperature 
scaling, mean temperature and seasonality) might be systematically 
biased and therefore their magnitudes are taken to be invariant with 
respect to the interval of averaging. On the basis of the 2e uncertain- 
ties for summer and winter temperature (Fig. 2a,c), we assessed the 
significance of dominant trends using Monte Carlo analysis (Fig. 2b,d; 
section on Trend analysis below). 


Seasonal preservation bias uncertainty. Unequal seasonal distribu- 
tion of snowfall could result in different magnitudes of diffusion for 
winter and summer amplitudes*. The seasonal temperature cycle also 
affects the magnitude of diffusion for all seasons. We used the Com- 
munity Firn Model (CFM)**”, a firn-evolution model with coupled firn 
temperature, firn densification and water-isotope modules, to test how 
seasonally weighted accumulation affects the diffusion of specified, 
hypothetical isotope records progressing from surface snow (8D,,,.y,), to 
consolidated snowpack in the firn (6D,,,,), to solid ice beneath the pore 
close-off depth (6D;.,). We applied the back-diffusion calculation (sec- 
tion on Diffusion corrections) to 8D; to estimate the original 5D,,oy. 
We then assessed how reconstructions of 6D,,,, could be misinter- 
preted as a result of different seasonal-accumulation weightings 
(Extended Data Fig. 2a,b). 

We performed five CFM runs using monthly time steps for accumula- 
tion, temperatureandisotopes (Extended DataTable1). Theseasonalcycle 
for 5D,,,ow is based on the mean amplitude in Fig. 1b (15.43%o). Five WAIS 
accumulation scenarios were tested on the basis of monthly accumula- 
tion from the regional climate model MAR3.6 (Modéle Atmosphérique 
Régional; ERA-Interim forced)", which spans the period January 1979 
to December 2017. The mean accumulation over the entire 39-yr 
period is 0.24745 m ice equivalent yr“, with about 1.6x as much snow in 
winter (April to September) as in summer (October to March). The five 
scenarios are as follows: (1) ‘constant’: identical accumulation for all 
months (0.0206 m ice equivalent month”; one-twelfth of the annual 
mean); (2) ‘cycle’: monthly accumulation equal to MAR monthly means; 
(3) ‘noise’: using the ‘cycle’ time series, we add noise to each time step in 
the ‘cycle’ series in the form ofa normal random variable of zero mean 
and the standard deviation for the month from MAR; (4) ‘random’: for 
each month, the accumulation is anormal random variable with mean 
and standard deviation equal to MAR monthly values; and (5) ‘loop’: the 
entire 39-yr MAR accumulation time series is repeated over and over 
again. For the temperature boundary condition, we used the mean tem- 
perature at a height of 2 m for 1979-2017 for each month predicted by 
MARto create an annual temperature cycle. We repeated this 12-month 
time series for the duration of the model runs. This method ensures 
that model runs, which are designed to test accumulation seasonal- 
ity, are not affected by interannual temperature variability, while also 
providing an estimate ofthe annual temperature cycle, which affects 
therate of isotope diffusion in the upper firn. 

Extended Data Fig. 2a,b shows the results for the ‘constant’ and ‘cycle’ 
cases. The diffusion-correction technique accurately reconstructs 
6D,,,, for summer and winter in the ‘constant’ snowfall scenario but 
underestimates summer values in the MAR ‘cycle’ scenario by about 
2.6%o, which is 11% of the full range of the observed WDC summer 
water-isotope values. Winter values are overestimated by only about 
0.6%, about 3% of the full winter range, since winter has 1.6x as much 
snowas summer. The ‘noise’, ‘random’ and ‘loop’ runs produce results 
within 0.3%o. These CFM experiments demonstrate that centennial 
trends in the summer and winter water isotopes of the order of a few 
per mil (%o) could arise from large changes in seasonal-accumulation 
weighting, whereas multimillennial trends >2.6%o are unlikely to be 
caused by seasonal accumulation and can therefore be interpreted as 
climate signals ofa different origin. For 1,000-yr averaging (as in Fig. 2), 
HadCM3 indicates seasonal-accumulation weighting (winter:summer) 


of 1.3 to 1.7 throughout the Holocene (Extended Data Fig. 2f), which 
yields a 1o uncertainty of 0.27%. based on the CFM testing criteria. 

To determine observationally if seasonal snowfall changed across 
the Holocene, we used measured black carbon (BC) concentrations, 
the only age-scale-independent impurity. BC data are available from 
0-2.5 ka and 6-11 ka (ref. ??). Seasonal-fire regimes in South America 
dominate BC concentrations at WDC, causing BC maxima and minima 
in autumn and spring, respectively? (Extended Data Fig. 2c). We split 
each year into two parts, characterized by rising or falling BC: BC, and 
BC, the depth intervals of rising and falling BC (Extended Data Fig. 2d). 
The duration of BC, is longer than BC, owing to source characteristics”, 
thus BC,/BC, < 1 (Extended Data Fig. 2e). The BC,/BC, ratio can change 
with time because of variability at the source, changes in atmospheric 
transport or seasonality of snow deposition. We observe little change in 
BC,/BC, resembling the multimillennial trends seen in WDC summers 
and winters (Extended Data Fig. 2e). Unless there are competing 
and exactly compensating effects in seasonality (the source change 
exactly cancels the depositional and transport change or other unlikely 
scenarios), the BC data provide evidence that changes in WDC seasonal 
snowfall were not large enough to affect our multimillennial climate 
interpretations. 


Intermittency of precipitation uncertainty. The episodic nature of 
snowfall creates an incomplete record of local climate variations, 
preventing interpretation of trends over short time intervals. We want 
to interpret isotopic variations averaged over a sufficiently long time- 
scale so that, to within a specified tolerance, trends are not likely to be 
random noise arising from the spread of distributions preserved inthe 
ice. Using distributions of reconstructed annual amplitudes (Fig. Ib and 
Extended Data Fig. 2g) for 1,000-yr windows throughout the Holocene, 
we conducted Monte Carlo resampling simulations to determine that 
250-yr averaging-lengths are needed to achieve a standard error of 1%o, 
corresponding to a mean amplitude-to-noise ratio of 15. For the time 
period with greatest variability, centred on 4 ka, the standard error for 
a1,000-yr average (as used in Fig. 2) is 0.52%o (Extended Data Fig. 2h). 
Because this is an amplitude uncertainty (rather than uncertainty 
associated with a season), we specify the lo uncertainties for summer 
and winter as half of 0.52%o. 


Isotope-temperature scaling and associated uncertainty. The con- 
version of isotopic values (1,000-yr averages) to temperature yields 
three curves for summer and three for winter: T omina 7,,, and To. Each 
curve is normalized to the value at 1 ka (as done in Fig. 2). The difference 
inthe 7,,,and 7_,,curves gives the louncertainty range *O; tO —Orscaler 
which are then added in quadrature to the ref. ? mean-temperature 
uncertainties, yielding the final uncertainty estimates shown in Fig. 2. 


Relationship between the annual mean and individual seasons. 
Using 1,000-yr and 300-yr averages of summer, winter and mean tem- 
perature (Extended Data Fig. 3c-f), we determined R^ values for summer 
and winter versus the mean. We then subtracted the 1,000-yr averages 
from the 300-yr averages to obtain residuals and then determined R? 
values again for summer and winter versus the mean (Extended Data 
Table 2). From 11to 0 ka, the high summer correlations for the 1,000-yr 
comparison indicate a strong association of the annual mean tempera- 
ture with the summer temperature at orbital timescales. At suborbital 
scales (300-1,000-yr residuals), neither the summer nor the winter 
alone explain much of the mean annual variability and the annual mean 
isarandom composite of the two seasons. If only 11-7 kais considered, 
winter variability explains more of the mean at submillennial scales. 


Trend analysis. To assess the significance of dominant trends in our 
reconstructed seasonal temperatures, we conducted a Monte Carlo 
analysis founded on the assumption that all possibilities for the 
unknown time-dependence of errors are equally likely. The essential 
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motivation for this approach is that we have determined the magni- 
tudes of uncertainties as a function of age but that we have no informa- 
tion about whether the errors in our reconstruction persist at similar 
values for long periods of time (exhibit a bias) or whether they fluctuate 
at high-frequency. 

We randomly generated a large number of alternative seasonal tem- 
perature histories governed by the uncertainties on 1,000-yr averages 
(Extended DataFig.3a,b), calculated the temperature trends for each 
alternative history over desired time intervals (such as 11-4 ka and 
4-0 ka) and compiled the results into frequency distributions from 
which probabilities can be calculated (Fig. 2b,d). Specifically, each 
alternative history deviates from the summer and winter temperature 
reconstruction by an amount that smoothly varies over time between 
random nodes whose values are a Gaussian random variable of zero 
mean and standard deviation for 1,000-yr averages at the age of the 
node. Thenumber of nodes and the age of each node are random vari- 
ables, uniformly distributed between 1and 11 nodes and 0-11 ka, respec- 
tively. A small number of nodes produce an alternative temperature 
history for which the bias is serially correlated for millennia, whereas 
alarge number of nodes producea history for which the bias is uncor- 
related from millennium to millennium. 


Seasonal moist energy balance model 

We used a simple global, zonal mean MEBM to calculate surface tem- 
peratures (Extended Data Fig. 4), accounting for TOA insolation, 
temperature-dependent longwave emission to space, temperature- 
dependent albedo to simulate brightening by snow and ice and hori- 
zontal atmospheric-heat transport treated as diffusion of near-surface 
moist static energy” ©. The model has a 2° spatial resolution, a sin- 
gle surface and single atmospheric layer and a surface-heat capacity 
based on the relative fraction of land and ocean surface in the zonal 
mean. Heat exchange between the surface and atmosphere layers 
arises from differences in blackbody radiation from each layer and 
sensible and latent heat exchanges proportional to temperature and 
specific humidity contrasts (assuming a constant relative humidity of 
80%), following bulk aerodynamic formulae. We calculated the annual 
TOA insolation-cycle at every latitude in 500-yr time slices from 0 to 
11 ka. For each time slice the model is run at 2 hour time resolution for 
30 model-years to reach equilibrium. 

Extended Data Fig. 4 compares the temporal evolution of sum- 
mer maximum and winter minimum heat divergence by the atmosphere 
(V-F) at the WDC site to the summer maximum and winter minimum 
site temperature and insolation. Although Holocene changes in V-F 
atthe WDC site correlate with insolation forcing, the magnitude of 
changes in maximum direct insolation are much larger than those in 
atmospheric-heat divergence. Further, the Holocene changes in sum- 
mertime heat divergence are ofthe wrong sign to cause net heating at 
the WDC site (positive divergence is an export of heat by the atmos- 
phere from the site). Heat transport in the Antarctic is convergent in 
the annual mean but divergent in mid-summer, as intense incoming 
insolation exceeds longwave emissions from the cold surface. 


HadCM3 model simulations 

Model setup. We used the fully coupled ocean-atmosphere model 
HadCM3 (refs. 9^5, v/HadCM3BM2.1, which well simulates tropical 
Pacific climate and its response to glacial forcing$6. Our simulations are 
snapshots at 1 kyr intervals over the last 11 ka (ref. ?), with time-specific 
boundary conditions of orbital forcing’, GHG concentration”, 
ice-sheet topography and sea level*^*7"7* We used three simulations: 
(1) only orbital forcing changes (ORBIT), with all other boundary con- 
ditions set to the pre-industrial; (2) orbital/GHG forcing with GLACID 
ice-sheet elevation history; and (3) orbital/GHG forcing with ICE-6G 
ice-sheet elevation history. Elevation histories are shown in Fig. 4a. 
Snapshot simulations were run for at least 500 yr with analysis made 
on the final 100 yr. Further snapshot simulations for 10 ka allowed 


us to decompose the role of different forcings, described in the fol- 
lowing sections. The large difference in forcings between 10 ka and 
the pre-industrial late Holocene epoch make this comparison most 
instructive. 


Summer climate. We examined the zonal mean at 80° S in simula- 
tions for hypothetical 10 ka worlds, by changing the boundary condi- 
tions to compare to the pre-industrial/late Holocene. These simu- 
lations are ‘10 ka ORBIT-only’; two runs with only ice sheets at 10 ka 
and pre-industrial settings otherwise, called '10 ka GLACID-only' and 
‘10 ka ICE-6G-only’; and two runs with all 10 ka forcings, called ‘10 ka 
GLACID-all' and ‘10 ka ICE-6G-all’. In ‘10 ka ORBIT-only’, reduced TOA 
shortwave radiation causes a large reduction in shortwave radiation at 
the surface (SW,) and consequent cooling. Downward longwave radia- 
tion (LW,) also decreases, probably because of atmospheric cooling. 
Sensible heat flux (SH,) to the surface is increased, indicating that the 
atmosphere and surface do not equally cool; one cause of this is the 
increased meridional heat convergence (-V-F). 

Changed ice sheets cause summertime cooling in both '10 ka 
GLACID-only' and ‘10 ka ICE-6G-only’, primarily via reduced LW,. 
Increased SW,, as a result of reduction of depth of the atmospheric 
column above the ice-sheet surface, counteracts the reduced LW, to 
someextent. (Reducingthe atmospheric column reduces SW absorp- 
tion and tends to cool the atmosphere, reducing LW,). Both ice-sheet 
scenarios also cause an increase in -V +F, partly counteracting sum- 
mertime cooling. 

Using all 10 ka forcings causes cooling through both LW, and SW,. 
The decrease in SW, is similar in 10 ka GLACID-all' and ‘10 ka ICE-6G-all’ 
and slightly smaller than in ‘10 ka ORBIT-only’, probably because the 
thinner atmospheric column reduces absorption. The decrease in 
LW, in 10 ka GLACID-all' and ‘10 ka ICE-6G-all' is larger than in ‘10 ka 
ORBIT-only’, ‘10 ka GLACID-only' and 10 ka ICE-6G-only’. This indicates 
the importance of feedbacks in the atmosphere. Heat convergence 
-V-F increases in both simulations indicating remote feedbacks, in 
addition to local feedbacks related to the amount of water vapour in 
the atmosphere. 

The preceding description of changes in the zonal mean in the 10 ka 
simulations compared to pre-industrial period holds for the entire 
Holocene epoch. Orbital forcing alone reduces SW, and LW, by roughly 
the same magnitude. With full forcing (including ice sheets), the reduc- 
tioninLW,isroughly three times the reduction in SW,. Considering an 
energy budget over the WDC site (79.467? S,112.085? W), mechanisms 
arethe same as for the zonal mean. Magnitudes of forcings change 
but reduced SW, still cools the surface, amplified by an LW, feedback 
dependent on ice-sheet size. 


Winter climate. During winter, SW, is no longer a factor as the sun 
is below the horizon, yet there is still surface warming caused by an 
increase in LW,. An increase in -V-F in ‘10 ka ORBIT-only’ warms the 
atmospheric temperature, increasing LW, and SHa. With an ice sheet 
imposed, the surface temperature cools. In both ‘10 ka GLACID-only' 
and ‘10 ka ICE-6G-only' there is a reduction in - V-F, reducing LW, and 
SH,. Whenall10 ka forcings are introduced, the change in temperature 
is smaller than for ice sheet-only runs. In ‘10 ka GLACID-all' we found 
no change in -V-F, LW, or surface temperature. This suggests that the 
increase in -V -F from orbital forcing is almost perfectly balanced by 
the change in -V-F from the ice-sheet configuration. 

The processes controlling heat transport over Antarctica are com- 
plicated and HadCM3 may not be able to simulate them perfectly. Our 
simulations indicated that remote processes during winter alter the 
heattransport, affecting atmospheric and surfacetemperatures. Rais- 
ing the topography of Antarctica tends to reduce such heat transport 
(Extended Data Fig. 5), producing an additional cooling on top of a 
pure lapse-rate effect”. This cannot, however, explain the prominent 
millennial-scale changes at about 9.2 and about 7.9 ka (Figs. 1e and 2c). 


Theintricacies ofinterpreting the early Holocene winter variability in 
West Antarctica necessitates further study. 


Sea ice. Sea ice changes may alter local energy fluxes from the ocean 
to the atmosphere. In HadCM3, sea ice extent changes across the 
Holocene. We used two analyses (Extended Data Fig. 6) to showthat sea 
ice is nota primary control on the surface temperature at WDC (80? S): 
(1) correlation analysis of sea ice changes and surface temperature and 
(2)atmosphere-only model simulations in which we specified individual 
changes in the model boundary conditions (including sea ice). 

We computed the dominant spatial patterns of sea ice variability 
using empirical orthogonal functions across all of the HadCM3 simu- 
lations (ALL) and individually for three subset simulations (ORBIT, 
GLACID and ICE-6G), from O to 11 ka. We projected the model-simulated 
seaicefor each individual time-slice simulation onto these patterns to 
compute the amplitude of sea ice variability in each simulation. The 
amplitude was compared to temperature at 80? S to understand how 
large-scale changes inthe seaice affect temperaturefor the months of 
December (summer) and July (winter) (Extended Data Fig. 6). 

In winter, we found negligible correlations between sea ice change 
and temperature in all sets of simulations (ALL: 0.02; ORBIT: 0.04; 
GLACID: 0.16; ICE-6G: 0.06). This suggests that winter sea ice is not 
animportant factor in determining the temperature at 80? S. In sum- 
mer, only the ORBIT simulation has meaningful correlations between 
temperature and sea ice variability (ALL: 0.59; ORBIT: 0.84; GLACID: 
—0.60; ICE-6G: 0.36). The sign of the correlation changes between 
simulations despite the sea ice change pattern being the same in all 
simulations. From this we concluded that sea ice is not a dominant con- 
trol on temperature at 80? S in summer. The correlation in the ORBIT 
simulations suggests that there may be some relationship between 
sea ice and temperature; we investigated this with atmosphere-only 
simulations. 

In atmosphere-only simulations at O ka and 10 ka, we specified the 
top ofthe atmosphere insolation, sea surface temperature (SST) and 
seaicefromthe ORBIT simulations. Over land areas and sea ice regions, 
the model calculates the surface temperature using the land-surface 
scheme inthe model. The atmosphere model is identical to the model 
used within the coupled model. We ran a series of experiments vary- 
ing the orbital configuration, SST or sea ice (summarized in Extended 
Data Table 3). 

The zonal mean of the change in the sea ice that we prescribed can be 
seen in Extended Data Fig. 6e and the change in the SST can be inferred 
from Extended Data Fig. 6f-h. Extended Data Fig. 6f shows that the 
atmosphere-only model replicates the change in temperature of the 
coupled model. Extended Data Fig. 6g shows that the effect of the 10 ka 
orbital configuration (‘Atmos_10k_insol’) is to cool Antarctica consider- 
ably by about 0.5 °C. North of 65° S there is no change in the surface 
temperature, primarily a response to the imposed SST and sea ice, 
which are the same in the ‘Control’ and Atmos 1Ok insol. Imposing 
the SST and sea ice from10 ka (‘Atmos_10k_ice SST’), we find very little 
change inthe surface temperature over Antarctica but there are some 
large changes in the surface temperature north of 70° S. Extended Data 
Fig. 6h shows the result of imposing 10 ka SST or 10 ka sea ice. The 10 ka 
SST (‘Atmos_10k_ SST’) tends to warm Antarctica, consistent with the 
large increases in SST north of 65? S. Changing sea ice (‘Atmos _10k ice’) 
tends to cool Antarctica. Both effects are small, approximately 0.1 °C 
and of opposite sign. This explains the small net change in the surface- 
temperature change over Antarctica when SST and sea ice are changed 
simultaneously, as shown by ‘Atmos_10k_SST_ice’. It should be noted 
that inthe coupled system a change in sea ice cannot be decoupled from 
a change in the SST, so not only is the effect of sea ice on the climate 
small, itis also probably associated with a compensating change in SST. 
From these simulations we concluded that sea ice and SST changes 
are not a dominant driver of the change in the surface temperature 
over Antarctica. 


Extended Data Fig. 6e shows that the change in sea ice at 10 ka in 
ORBIT is much larger than the change in either GLACID or ICE-6G. The 
ORBIT simulations do not account for all ofthe changes in the boundary 
conditions at 10 ka andare therefore less realistic than either ICE-6G or 
GLACID. Because ICE-6G and GLACID both show much smaller changes 
inthe sea ice and SST than ORBIT, we expect that in reality there is also 
amuch smaller change in the sea ice and SST than in ORBIT. We thus 
concluded that sea ice has a small impact on the temperature at 80? S 
in summer. 

We also performed a similar analysis ofthe winter season (not shown). 
Wefoundthatthe atmosphere-only model does not compare well with 
thecoupled model, simulating very little change in the surface tempera- 
ture. Doing aterm-by-term decomposition of the atmosphere model 
isnot, therefore, particularly useful asittells us more about the model 
rather than the physical climate. The failure of the atmosphere-only 
modelto capture the changes at 10 kasuggests that the importance of 
the SST and sea ice is in their day-to-day coupling with the atmosphere 
and not in any long-term mean change in this season. 


Estimating elevation changes 

Temperatures simulated by HadCM3 for ORBIT provide a control sce- 
nario against which observations can be compared to identify the signal 
of elevation change. For a chosen time interval, the net reconstructed 
warming AT, exceeds that of ORBIT by an amount A7; - AT,. This can 
be compared to the effective lapse rate (AT - AT9)/AZ,, defined by a 
HadCM3 simulation including topographic change (GLACID or ICE-6G 
models) and all forcings, for model warming AT,, and model eleva- 
tion decrease AZ,. Specifically, the estimated elevation decrease is 
AZ, = AZy((ATa 7 ATo)/(ATy - ATo)). Summer and winter reconstruc- 
tions offer two separate assessments, for which we calculate the alge- 
braic average. 

Accounting for uncertainties in AT; requires recognizing that uncer- 
tainties of summer and winter reconstructions are not independent, 
whilealso recognizing that they emerge from two independent sources: 
uncertainty in the mean annual temperature history (calculated in ref. ??) 
and uncertainty in the seasonal amplitude (calculated in the present 
study). In general, the uncertainty of seasonal temperature at a speci- 
fiedtimeisthe quadrature sum of annual and amplitude uncertainties, 
which gives Ic and 2c uncertainties for a given season and time. How- 
ever, if the true value of annual temperature is shifted by an amount ao 
from the nominal reconstruction, this must be true for both summer 
and winter. And if the true value of amplitude is shifted by an amount 
Bo from the nominal reconstruction, the temperature shift must be 
+Boin one season but -$o in the other. 

To define bounding cases on elevation change in a specified time 
interval, we calculated the maximal (or minimal) temperature change 
AT, for a season by differencing the upper (or lower) limit at one end 
of the interval with the lower (or upper) limit at the other end and also 
calculating the corresponding A7, for the opposite season required by 
the correlated errors. The elevation decreases AZ, were then calculated 
by comparison to HadCM3 simulations, as specified previously and 
the summer and winter values averaged. This process was completed 
four times for each time interval and HadCM3 model, corresponding 
to four different initial AT, (maximum and minimum AT} for summer 
and maximum and minimum AT, for winter) and the most extreme case 
taken as the result (this proved to be the one starting with maximum 
summer AT,). Table 1 lists results for two time intervals and the two 
HadCM3 simulations with variable topography. 


Data availability 

The WDC water-isotope datasets analysed during the current study are 
available inthe US Antarctic Program Data Center (USAP-DC) repository, 
https://doi.org/10.15784/601274 and https://doi.org/10.15784/601326. 
The impurity datasets analysed during the current study are available 
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in USAP-DC repository, https://doi.org/10.15784/601008. The data 
generated inthis study areavailable in the USAP-DC repository, https:// 
doi.org/10.15784/601603, including raw and diffusion-corrected water 
isotopes, seasonal water isotopes (maximum summer and minimum 
winter values) and seasonal temperature reconstructions. Source data 
are provided with this paper. 


Code availability 


The MATLAB code used for the diffusion correction of water-isotope 
data and the subsequent selection of seasonal extrema (summer 
and winter) are available online at Zenodo, https://doi.org/10.5281/ 
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Extended Data Fig. 2 | Precipitation uncertainties. Uncertainties for 
seasonally weighted accumulation are shown in a-f, and for precipitation 
intermittency in g-h. a, The diffusion envelope of CFM output data (50-yr avg.), 
based onan input sine wave with f= lyr‘ and amplitude = 15.43%. The original 
(input) amplitude of the signal (black, dotted-dashed lines) decreases as time 
passes due to diffusion and downward advection of the firn, as shown by the 
decay of the maximum (red) and minimum (blue) lines, while the mean values 
of the ‘constant’ and ‘cycle’ (black solid and dashed lines, respectively) 


scenarios donot change and are dependent onthe seasonal weighting of snowfall. 


b, Diffused CFM output data from beneath the pore close-off depth (>200 yr) 
(black lines; smaller amplitude), with diffusion-corrected data shown with grey 
lines (larger amplitude). Red circles are the annual maximum value, and blue 
circles the annual minimum, selected using the same algorithm as Fig. 1a. 

c, Azoomofblack carbon (BC) concentrations at -6.5 ka??5?, The maxima 


(red circles) and minima (blue circles) can be used to separate approximate 
depth intervals corresponding to winter (BC,) and summer (BC,); vertical blue 
lines correspond to nominal January 1, as defined by the peak of nssS/Na”. 

d, The140-yr averages for BC, (blue) and BC, (red). The grey lineis WDC annual 
accumulation”, orange circles are BC, + BC,, which should equal annual 
accumulation. e, Black carbon seasonality BC;/BC, (black), based on (d). 

f, Accumulation seasonality for HadCM3 seasonal snowfall (red line) compared 
to the range of seasonality tested using the CFM (dashed blue lines) and 
modern MAR seasonality® (blue diamond). g, Distribution of annual amplitudes 
for water isotopes for a1,000-year window centred at 4 ka. h, Standard errors 
are determined for1,000-realizations of random sampling of the distribution 
in (g) to determine a standard deviation of the residuals of the true mean minus 
the mean ofthe random sampling. 
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MEBM annual results for the mean, maximum, and minimum are shown in d-l. to heat convergence at the site. 
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Positive values all indicate a surface or atmospheric warming. Variables include 
surface temperature (T, in °K), latent heat (LH in Wm), sensible heat to the 
surface (SH,in Wm?), shortwave radiation (SW in Wm?), downward LW 
radiation (LW,in Wm”), and change in heat transport (V-F in 107 W). 
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average. e, Change in sea ice fraction for coupled-model simulations from 

10 katoO ka. f, Change in surface temperature between 10 ka and O ka for the 
coupled-model simulations (dotted line) and atmosphere-only simulations 
(solid line). g, Change in surface temperature for atmosphere-only runs from 
10 katoO ka. h, Change in surface temperature for atmosphere-only runs 
from10 kato 0 ka. See Extended Data Table 3 for descriptions of the model 
simulations used in panels (g) and (h). 
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1 kyr resolution, HadCM3 ORBIT vs. WDC Temperature 
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MEBM Winter, R 2. 0.05, p-value 0.30 


eo 
p 


o 


MEBM Winter Temperature Anomaly (°C) 
o 
N 


0 0.2 0.4 0.6 0.8 
WDC Winter Temperature Anomaly (°C) 


500 yr resolution, MEBM vs. WDC Temperature 
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Extended Data Fig. 7 | Model results vs. WDC temperatures. Coefficient of determination and p-values for comparison of HadCM3 (1-kyr resolution, n =12) or 
MEBM (0.5-kyr resolution, n = 23) model results with WDC summer and winter temperatures. 


Extended Data Table 1| CFM simulation inputs 


Month | b (m ice eq. yr?) &(miceeq.yr)| TCO SD (%o) 
1 0.139 0.088 -17.5 -250.22 
2 0.169 0.072 -24.9 -252.48 
3 0.267 0.010 -31.5 -258.64 
4 0.293 0.121 -34.6 -267.05 
5 0.345 0.137 -34.9 -275.46 
6 0.311 0.136 -362 -281.61 
7 0.310 0.111 
8 0.298 0.116 
9 0275 0.112 
10 0.271 0.136 
1 0.164 0.079 
12 0.128 0.071 


Meanb and standard deviation Og of each month's accumulation rate and monthly mean temperature at WAIS Divide for 1979 to 2017 predicted by MAR? and the isotope values used for each 
month during the CFM simulations. 
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Extended Data Table 2 | Seasonal vs. mean correlations 


1,000-year 300-year 300 to 1,000- 
11-0 ka 1 
averages averages year residuals 


mean 


mean 


1,000-year 300-year 300 to 1,000- 
11-7 ka k 

averages averages year residuals 
dd 
mean 
uu 
mean 


R? values for 1,000-year and 300-year averages of summer and winter vs. mean, as well as for residuals. 


Extended Data Table 3 | Atmosphere-only climate model experiments 


SST Seaice | TOA insolation 


Control Oka Oka Oka 
Atmos 10k 10ka 10ka 10ka 
Atmos_10k_ice/SST 10ka 10ka Oka 
Atmos_10k_insol Oka Oka 10ka 
Atmos_10k_ice Oka 10ka Oka 
Atmos 10k SST 10k Oka Oka 


Variations of boundary conditions in each experiment. 


